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ABSTRACT 

We present updates to PRISM, a photometric transit-starspot model, and GEMC, a hybrid op¬ 
timisation code combining MCMC and a genetic algorithm. We then present high-precision 
photometry of four transits in the WASP-6 planetary system, two of which contain a starspot 
anomaly. All four transits were modelled using PRISM and GEMC, and the physical properties 
of the system calculated. We find the mass and radius of the host star to be 0.836 zb 0.063 Mq 
and 0.864 zb 0.024 Rq, respectively. For the planet we find a mass of 0.485 zb 0.027 Mj^p, 
a radius of 1.230 zb 0.035 Rjup and a density of 0.244 zb 0.014 pjup- These values are con¬ 
sistent with those found in the literature. In the likely hypothesis that the two spot anomalies 
are caused by the same starspot or starspot complex, we measure the stars rotation period 
and velocity to be 23.80 zb 0.15d and 1.78 zb 0.20kms“^, respectively, at a co-latitude of 
75.8°. We find that the sky-projected angle between the stellar spin axis and the planetary 
orbital axis is A = 7.2° zb 3.7°, indicating axial alignment. Our results are consistent with and 
more precise than published spectroscopic measurements of the Rossiter-McLaughlin effect. 
These results suggest that WASP-6 b formed at a much greater distance from its host star and 
suffered orbital decay through tidal interactions with the protoplanetary disc. 

Key words: planetary systems — stars: fundamental parameters — stars: spots — stars: 
individual: WASP-6 — techniques: photometric 
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1 INTRODUCTION 

At present'^ a total of 1890 planets outside of our own solar sys¬ 
tem are listed in the authoritative catalogue of Schneider et al. 
(2011). Of these approximately two thirds have been discovered 
from ground-based (e.g. SuperWasp: Pollacco et al. 2006; HAT: 
Bakos et al. 2004) or space-based (CoRoT: Baglin et al. 2006; Ke¬ 
pler: Borucki et al. 2010) transit surveys, and later confirmed by use 
of the radial velocity technique (Butler et al. 1996, 1999; Queloz 
et al. 2000). Many more candidate exoplanets have been listed in 
the literature, mainly from the Kepler satellite survey which has 
also detected several Earth-size planets in the habitable zone (HZ) 
of their parent star, indicating new worlds with mass and size simi¬ 
lar to our own Earth (Borucki et al. 2012; Borucki et al. 2013). 

During a planetary transit, the planet follows a path (called 
the transit chord) across the surface of the stellar disc and can be 
used to probe changes in brightness on the stellar surface (Silva 
2003). Starspots have different temperatures to the surrounding 
photosphere, so emit a different amount of flux. Because photom¬ 
etry measures the change in intensity as a function of time, the oc- 
cultation of a starspot by the planet causes an anomaly in the light 
curve (Silva 2003). The anomaly is either an increase or decrease 
in the amount of light received from the star. If the starspot is a cool 
spot then the amount of light will increase when the planet crosses 
the starspot (Rabus et al. 2009; Pont et al. 2007; Winn et al. 2010b). 
If the starspot is a hot spot (e.g. a facula) then the amount of light 
will reduce when the planet occults the spot. 

At present, when a light curve of a transiting exoplanet is 
observed to have a starspot anomaly, the transit and the spot are 
generally modelled separately (e.g Desert et al. 2011; Maciejew- 
ski et al. 2011; Nutzman et al. 2011; Sanchis-Ojeda et al. 2011). 
First, a transit model is fltted to the datapoints not affected by the 
starspot anomaly. Then the spot-affected residuals versus the best¬ 
fitting model are modelled using a Gaussian function (e.g. Sanchis- 
Ojeda et al. 2011; Sanchis-Ojeda & Winn 2011). This method ne¬ 
glects the fact that the starspot affects the entire transit shape and 
not just the section where the planet crosses the spot (Ballerini et al. 
2012). Carter et al. (2011) use the idea that a starspot on the stellar 
disc will affect the transit depth to explain the observed changes 
in transit depth for GJ 1214. This is due to the change in the star’s 
brightness in its long-term light curve due to starspots rotating on 
and off the stellar disc. 

The transit depth is not the only property of a transit light 
curve that the starspot affects: it also affects the determination of 
the measured stellar mean density, stellar radius, orbital inclination 
and limb darkening (LD) coefficients (Ballerini et al. 2012). The 
LD coefficients depend on wavelength: because a starspot has a dif¬ 
ferent temperature compared to the surrounding photosphere, it has 
a different spectral energy distribution and thus different LD coeffi¬ 
cients. Therefore the application of a LD law with a single set of co¬ 
efficients to the entire stellar surface causes a bias in the modelling 
process (Ballerini et al. 2012). The difference in LD coefficients 
between the spot and the photosphere can be as much as 30% in 
the UV. The effects on the measured stellar radius and orbital incli¬ 
nation of the system are artifacts from errors in the measured plane¬ 
tary radius, which is derived from the transit depth. A change in the 
measured planetary radius must be compensated for by a change in 
the measured stellar radius or semimajor axis in order to retain the 
same transit duration. Starspots can also affect the measured transit 
midpoint (Sanchis-Ojeda et al. 2011; Barros et al. 2013) and create 

^ (http://exoplanet.eu) accessed on 2015/02/20 


false positives in transit timing measurements. Sanchis-Ojeda et al. 
(2011) calculated that a starspot anomaly in a transit of WASP-4 
with an amplitude of 0.3 to 0.5 mmag could produce a timing noise 
of five to ten seconds. 


1.1 Introducing WASP-6 

The transiting planetary system WASP-6 was discovered by Gillon 
et al. (2009) using photometry from the WASP-South telescope. 
They determined an orbital period of P = 3.361 days for the 
planet WASP-6 b. Dedicated photometric observations were then 
performed in the i' band using the 2-m Faulkes Telescope South 
(FTS) and in a broad V+R band using the RISE instrument (Steele 
et al. 2008) on the 2-m Liverpool Telescope (LT). 

Radial velocity (RV) measurements were obtained using two 
spectrographs: CORALIE on the L2-m Euler telescope (Baranne 
et al. 1996; Queloz et al. 2000) and HARPS on the ESO 3.6- 
m telescope (Mayor et al. 2003). Gillon et al. (2009) determined 
the stellar mass and radius to be M* = O.SS^Q Qg M© and 
P* = 0.870j)Q Q 3 g Rq, respectively. They found the planetary 
mass and radius to be Mp = O.SOSlQ Qgg Mjup and Rp = 
1.224lg Qg 2 Rjup- They also determined a value for the projected 
stellar rotational velocity of w sin 1 = 1.4 ± 1.0 kms“^ from mea¬ 
surements of line widths in the HARPS spectra with a macrotur¬ 
bulence (Umac) value of 2kms“^. They noted that if a value of 
Vmac = 0 km s”*^ is used then v sin 7 = 3.0 ± 0.5 km s“^, while if 
Vmac became slightly larger than 2 km s“^ then v sin 7 would drop 
to zero. From their RVs Gillon et al. (2009) measured the Rossiter- 
McLaughlin (RM) effect. They found that the system is in align¬ 
ment with a sky-projected spin orbit alignment, A = lllis 

The spectral analysis of 11 WASP host stars by Doyle et al. 
(2013) included WASP-6 A. Doyle et al. (2013) derived new val¬ 
ues for the stellar mass and radius of A7* = 0.87 ± 0.06 M© and 
P* = 0.77 ± 0.07 Rq, in agreement with those of Gillon et al. 
(2009). Doyle et al. (2013) determined Vmac ~ 1.4 ± 0.3kms“^ 
and usin7 = 2.4 ± 0.5kms“^, and an effective temperature of 
Teff = 5375 ± 65 K. 

An optical transmission spectrum for WASP-6 has been con¬ 
structed using multi-object differential spectrophotometry with the 
IMACS spectrograph on the Magellan Baade telescope (Jordan 
et al. 2013). The observations comprised 91 spectra covering 480- 
860 nm. The analysis yielded a mostly featureless transmission 
spectrum with evidence of atmospheric hazes and condensates. 
Most recently, Nikolov et al. (2014) used the Hubble Space Tele¬ 
scope to perform transmission spectroscopy of WASP-6, and found 
a haze in the atmosphere of WASP-6 b. They also determined a ro¬ 
tational modulation of P-ot = 23.6 ± 0.5 d for WASP-6 A. 


2 UPDATES TO PRISM & GEMC 

A code written in IDL^ called PRISM (Planetary Retrospective Inte¬ 
grated Star-spot Model) was developed to model a starspot anomaly 
in transit light curves of WASP-19 (see Tregloan-Reed et al. 2013). 
PRISM uses a pixellation approach to represent the star and planet 
on a two-dimensional array in Cartesian coordinates. This makes it 
possible to model the transit, LD and starspots on the stellar disc 

^ The acronym IDL stands for Interactive Data Language and is a 
trademark of Exelis Visual Information Solutions. For further details see 
http://WWW.exelisvis.com/ProductsServices/IDL.aspx. 
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simultaneously. LD was implemented using the standard quadratic 
law. PRISM uses the ten parameters given in Table 1 to model the 
system, where the fractional stellar and planetary radii are defined 
as the absolute radii scaled by the semimajor axis (r*,p = i?*,p/a). 

A new optimisation algorithm called GEMC (Genetic Evolu¬ 
tion Markov Chain) was also created alongside PRISM to help im¬ 
prove the efficiency of finding a global solution in a rugged param¬ 
eter space compared to conventional MCMC algorithms (Tregloan- 
Reed et al. 2013). GEMC is a hybrid between an MCMC and a ge¬ 
netic algorithm® and is based on the Differential Evolution Markov 
Chain (DE-MC) put forward by Ter Braak (2006). During the 
‘burn-in’ stage GEMC runs N chains in parallel and for every gen¬ 
eration each chain is perturbed in a vector towards the current best¬ 
fitting chain. Once the burn-in stage has been completed GEMC 
switches to a conventional MCMC algorithm (each chain used in 
the burn-in begins independent MCMC runs) to determine the pa¬ 
rameter uncertainties (see Tregloan-Reed et al. 2013). 

Since the development of both PRISM and GEMC other authors 
have used the codes to not only help ascertain the photometric pa¬ 
rameters of a transiting system but to also derive the parameters 
of the starspots observed in transit light curves (e.g. Mancini et al. 
2013, 2014; Mohler-Fischer et al. 2013). Beky et al. (2014) used 
PRISM to help calibrate their semi-analytic transit-starspot model 
SPOTROD. 

Before using both PRISM and GEMC in modelling the WASP- 
6 system, it was decided to make a few improvements'^. The orig¬ 
inal version of PRISM assumed a circular orbit, as most transiting 
planets either have a circular orbit or lack a measurement of the or¬ 
bital eccentricity. However, Gillon et al. (2009) found that the orbit 
of WASP-6b has a small orbital eccentricity of e = O.OSdto'oJ®, 
with an argument of periastron cj = 97 .4^®3®2 degrees. As a con¬ 
sequence PRISM was extended to allow for eccentric orbits, e and 
u> have been set to roam within the physically bounded ranges of 
0 ^ e 1 and 0° a; 360°. A Gaussian prior is used to help 
constrain the parameter values close to the expected values found in 
the literature. The logic behind using a Gaussian prior stems from 
the fact that it is not possible to ascertain these values from pho¬ 
tometry alone (due to only observing a small fraction of the orbit) 
unless an occultation is observed (Kipping et al. 2012). Because 
we have the knowledge of where the values of e and cj should lie 
and that they have an effect on the other system parameters (in par¬ 
ticular i and r*), it is imperative to examine every potential solu¬ 
tion selected from a Gaussian probability distribution of e and uj 
to accurately estimate the uncertainties in all of the other system 
parameters. 

It was shown by both Silva-Valio et al. (2010) and Mohler- 
Fischer et al. (2013) that in some cases there can be more than one 
starspot anomaly in a single transit light curve. While PRISM was 
originally designed to model multiple starspots, the static coding of 
GEMC made it only possible to fit for either a single starspot or a 
spot free stellar surface. To facilitate further work GEMC was modi¬ 
fied to fit for multiple starspots. This was accomplished by allowing 
the initial reading of the input file to be dynamic, so GEMC can de¬ 
termine the number of starspots to be fitted based on the number 
of parameters used. This can be done by adding multiple spot pa¬ 
rameter ranges in the input file. It is possible to fix the position of 

® A genetic algorithm mimics biological processes by spawning successive 
generations of solutions based on breeding and mutation operators from the 
previous generation. 

^ The new versions of both PRISM and GEMC are available from 
http://www.as tro.keele.ac.uk/^Jtr 


a starspot and therefore assign starspots to sections of the stellar 
disc where they will not be occulted by the planet, thus allowing 
investigations of the effects of unocculted starspot on transit light 
curves. 

PRISM was designed to use the pixellation approach, and to 
maintain numerical resolution was hard-coded to set the planetary 
radius at 50 pixels. The host star’s radius in pixels was scaled ac¬ 
cordingly based on the input parameters. The new version now al¬ 
lows users to set the size of the planetary radius in pixels. This 
makes it possible to reduce the amount of time required to com¬ 
plete each model iteration, at the cost of numerical resolution (see 
Section 2.1 for more details). In tests using a planet radius set at 
15 pixels it took PRISM and GEMC approximately 13 s to model 
a single generation of 256 solutions using synthetic data, which 
equates to approximately 0.05 seconds per iteration. For compari¬ 
son, a planetary radius of 50 pixels results in approximately 0.47 s 
per iteration.® 

To increase the efficiency of determining the parameter uncer¬ 
tainties, the MCMC component of GEMC was replaced with DE- 
MC (Ter Braak 2006). DE-MC combines the genetic algorithm dif¬ 
ferential evolution (DE) (Price & Stom 1997; Stom & Price 1997) 
with MCMC. The combination of DE and MCMC is used to solve a 
problem in MCMC by determining the orientation and the scale of 
the step sizes. Adaptive directional sampling in MCMC does solve 
the orientation problem, but not the scale (Ter Braak 2006). DE- 
MC works by creating a population of MCMC chains whose start¬ 
ing points are initialised from overdispersed states and instead of 
letting the chains run independently and checking for convergence 
(e.g. Gelman & Rubin 1992) they are instead run in parallel and 
learn from each other. The perturbation steps taken by each chain 
are given by Eq. 1. Assuming a d-dimensional parameter space and 
using N chains then the population X is a x d matrix, with the 
chains labelled xi, X 2 ,... xjv. Therefore the proposal vector Xp is 
generated by: 

Xp = Xi -I- 7 (Xfli - Xfl 2 ) + e (1) 

where Xi is the current i*® chain, 7 is the scale factor calculated 
from 7 = 2.4/v/2d (Ter Braak 2006), xui and Xfl 2 are two ran¬ 
domly selected chains and e is drawn from a symmetric distribu¬ 
tion with a small variance compared to that of the target. Xp is then 
tested for fitness and if accepted it is used as the next step in Xi. 

After the ‘burn in’ stage of a MCMC chain, determining the 
required step size to allow a 20-25% acceptance rate can be diffi¬ 
cult. For a transit light curve altering the orbital inclination, i, by 
0.05% should only cause a small increase in but u 0.05% alter¬ 
ation in the transit midpoint. To, could cause a large increase in x^- 
DE-MC overcomes the problem with the scale of the step sizes by 
using the clustering of the chains around the global solution after 
the ‘burn in’: the difference vector between two randomly selected 
chains will contain the individual scale for each parameter (e.g. 
0.05% for i and 0.00001% for To). Ter Braak (2006) argues that 
DE-MC is a single Al-chain that is simply a single random walk 
Markov Chain in a x d dimensional space. 

The use of DE-MC in the exoplanet community is increas¬ 
ing, especially for models involving a large number of parameters. 
For example, models of transiting circumbinary planets can con¬ 
tain over 30 parameters (e.g. Doyle et al. 2011; Orosz et al. 2012; 
Welsh et al. 2012; Schwamb et al. 2013). To accurately estimate the 
parameter uncertainties the MCMC component of GEMC required 

® These tests were performed on a 2.4 GHz quad-core laptop. 
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Table 1. Original and recovered parameters from a simulated transit light curve using either 15 or 50 pixels for the planetary radius, plus the interval within 
which the best fit was searched for using GEMC. 


Parameter 

Symbol 

Original value 

Search interval 

Recovered value 
rp = 50 pixels 

Recovered value 

Tp = 15 pixels 

Radius ratio 

'■pA* 

0.15 

0.05 to 0.30 

0.1496 ±0.0013 

0.1498 ±0.0011 

Sum of fractional radii 

rs + rp 

0.25 

0.10 to 0.50 

0.2486 ± 0.0024 

0.2512 ±0.0026 

Linear LD coefficient 

Ui 

0.3 

0.0 to 1.0 

0.291 ±0.104 

0.281 ±0.114 

Quadratic LD coefficient 

U2 

0.2 

0.0 to 1.0 

0.192 ±0.042 

0.189 ±0.039 

Orbital Inclination (degrees) 

i 

85.0 

70.0 to 90.0 

85.16 ±0.46 

85.29 ± 0.44 

Transit epoch (Phase) 

To 

0.015 

-0.50 to 0.50 

0.01494 ± 0.00011 

0.01502 ± 0.00010 

Longitude of spot (degrees) 

e 

30.0 

-90.0 to +90.0 

30.50 ± 1.17 

30.47 ± 1.21 

Co-latitude of spot (degrees) 

<t> 

65.0 

0.0 to 90.0 

64.51 ±5.83 

64.17 ±5.55 

Spot angular radius (degrees) 

r'spot 

12.0 

0.0 to 30.0 

12.73 ± 2.00 

12.33 ± 1.87 

Spot contrast 

Pspot 

0.8 

0.0 to 1.0 

0.797 ± 0.057 

0.781 ± 0.061 


10® function iterations (10 chains each of 10® steps). The DE-MC 
component requires approximately 2 x 10® function iterations (128 
chains each of 1500 steps (e.g. Welsh et al. 2012)). This equates to 
a five-fold reduction in the amount of computing time required to 
fit a transit light curve. When using a set of synthetic transit data 
it took GEMC approximately 5.4 days to fit the data using a planet 
radius of 50 pixels coupled with the MCMC component. The use of 
a planet radius of 15 pixels combined with the DE-MC algorithm 
resulted in GEMC taking only 2.7 hours to fit the same data. 

2.1 Forward simulation of synthetic data 

The modifications to PRISM and GEMC were validated by mod¬ 
elling simulated transit data containing a starspot anomaly. For this 
test PRISM was used to create multiple simulated transits with a 
range of parameters. Noise was then added to the light curves so 
that the rms scatter between the original simulated light curves and 
the light curves with added noise was « 500 ppm. Other levels of 
noise where also used in similar tests. This was to approximate a 
realistic level of noise found in transit light curves observed using 
the defocused photometry technique. Error bars were then assigned 
to each datapoint to give the original noise-free model a reduced chi 
squared value of = 1- 

Once a simulated transit light curve had been created, GEMC 
and PRISM were used in an attempt to recover the initial input pa¬ 
rameters. Different values for the planetary pixel radius were also 
used to test for numerical resolution. Table 1 shows the results for 
one of the tests using both rp = 50 and Cp = 15 pixels, while Fig. 1 
shows the simulated transit light curve together with the original 
and recovered models for the same test using Pp = 50 pixels. 

From studying both Table 1 and Fig. 1, it can be seen that the 
recovered parameter values agree with the original values within 
their Icr uncertainties. Interestingly, the rms scatter of the recov¬ 
ered model was found to be 499 ppm while the rms scatter of the 
original model is 511 ppm. This showed that GEMC not only ex¬ 
plored the large parameter search space but also scanned the local 
area around the global solution to find the best possible fit® to the 
simulated data. This result is expected, and a testament to an opti¬ 
misation algorithm designed to find the lowest achievable xt (the 
recovered solution in this case had a xJ = 0.94) in a given param¬ 
eter space. Similar results were found on all the simulation tests 


° This best fit is in fact a phantom solution generated by the addition of 
noise. 



-0.10 -0.05 0.00 0.05 0.10 

Phase 


Figure 1. Recovered and original models to simulated transit data created 
by PRISM and recovered by GEMC and PRISM. The residuals are shown at 
the bottom. The model was calculated with rp = 50 pixels. 


and show that both GEMC and PRISM are capable of accurately and 
precisely determining the properties of transit light curves. 

The recovered parameter values from setting rp = 50 and 
Vp = 15 pixels also agree within their Icr uncertainties (see Ta¬ 
ble 1). The scale of the la uncertainties for when Vp = 15 are com¬ 
parable in scale to that of the la uncertainties for when rp — 50. 
This indicates that using a smaller number of pixels for the plane¬ 
tary radius (this reduction depends on the number of datapoints and 
the overall scale of the system being modelled) has little effect on 
the numerical resolution of the determined parameters or their as¬ 
sociated uncertainties. However, using a smaller number of pixels 
for the planetary radius does affect the smoothness of the plotted 
best-fit model. It is therefore advisable that, once the best-fitting 
parameters have been found, GEMC is used again with the param¬ 
eters fixed at the best-fitting values and with rp = 50 to calculate 
a smooth best-fitting model. These tests showed that it is possible 
to obtain precise results and correctly estimated parameter uncer¬ 
tainties, whilst, using a planetary pixel radius of less than 50. There 
are, though, some values which should not be used. For example 
in tests using rp = 5 the parameter uncertainties were heavily un- 
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derestimated, due to numerical noise in the model. By making the 
planet only 10 pixels across, the numerical resolution decreases to 
the point where adverse effects can be seen in the results and un¬ 
certainties. 


3 OBSERVATIONS AND DATA REDUCTION 

Four transits of WASP-6 were observed on 2009/06/26, 
2009/08/02, 2009/08/29 and 2010/07/31 by the MiNDSTEp con¬ 
sortium (Dominik et al. 2010) using the Danish 1.54 m telescope 
at ESO’s La Silla observatory in Chile. The instrument used was 
the DFOSC imager, operated with a Bessell R filter. In this setup 
the CCD covers a field of view of (13.7')^ with a pixel scale 
of 0.39” pixel“^. The images were unbinned but windowed for 
faster readout, resulting in a dead time between consecutive im¬ 
ages of between 22 and 35 s. The exposure times were 80-120 s. 
The Moon’s brightness and distance to the target star is given in 
Table 2. The telescope was defocused and autoguiding was main¬ 
tained through all observations. The amount of defocus applied 
caused the resulting PSFs to have a diameter of 86 pixels for the 
night of 2009/06/26, 32 pixels for the night of 2009/08/02, 44 
pixels for the night of 2009/08/29 and 37 pixels for the night of 
2010/07/31. 

We reduced the data in an identical fashion to Southworth 
et al. (2009a,b). In short, aperture photometry was performed with 
an IDL implementation of DAOPHOT (Stetson 1987), and the aper¬ 
ture sizes were adjusted to obtain the best results (see Table 2). A 
first order polynomial was then fitted to the outside-transit data 
whilst simultaneously optimising the weights of the comparison 
stars. The resulting data have scatters ranging from 0.591 to 1.215 
mmag per point versus a transit fit using PRISM. The timestamps 
from the fits files were converted to BJD/TDB. An observing log is 
given in Table 2 and the final light curves are plotted in Fig.2. 


4 DATA ANALYSIS 

All four transits were modelled using PRISM and GEMC. To do this 
a large parameter search space was selected to allow the global 
best fit solution to be found. As discussed in Tregloan-Reed et al. 
(2013), the ability of GEMC to find the global minimum in a short 
amount of computing time meant that it was possible to search a 
large area of parameter space to avoid the possibility of missing 
the best solution. The parameter search ranges used in analysing 
the WASP-6 datasets are given in Table 3. We modelled the two 
datasets containing a starspot anomaly independently, in order to 
obtain two sets of starspot parameters. This helps the investigation 
of whether the two anomalies are due to the same starspot (see Sec¬ 
tion 5). 

The separate models of the four datasets of WASP-6 have pa¬ 
rameters which are within la of each other (Table 3). Ballerini et al. 
(2012) noted that starspots can affect the LD coefficients by up to 
10% in the R band. This is not seen in the WASP-6 data, unlike 
in the transit data of WASP-19 (Tregloan-Reed et al. 2013). The 
scatter around the weighted mean is yf, = 0.149 for the linear co¬ 
efficient and 0.355 for the quadratic coefficient. The error bars on 
the LD coefficients are too large to allow the effects of starspots to 
be detected. This is due to the lower quality of the data compared to 
WASP-19 (Tregloan-Reed et al. 2013). The combined best-fit LD 
coefficients are also in agreement within their la uncertainties with 
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Figure 2. The four light curves of WASP-6 presented in this work, in the 
order presented in Table 2. Times are given relative to the midpoint of each 
transit. 


the theoretically predicted values for WASP-6 A of ui = 0.4125 
and U 2 = 0.2773 (Claret 2000). 


4.1 Photometric results 

The final photometric parameters for the WASP-6 system are given 
in Table 4 and are weighted means together with their la uncertain¬ 
ties of the results from the four individual fits. Fig. 3 compares the 
light curves to the best-fitting models, including the residuals. 

The available times of mid-transit for WASP-6 were collected 
from the literature (Gillon et al. 2009; Dragomir et al. 2011; Sada 
et al. 2012; Jordan et al. 2013; Nikolov et al. 2014). All timings 
were converted to the BJD/TDB timescale and used to obtain an 
improved orbital ephemeris: 

To = BJD/TDB 2 454 425.02180(11) -f 3.36100208(31) x E 
where E represents the cycle count with respect to the reference 
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Table 2. Log of the observations presented for WASP-6. is the number of observations. ‘Moon ilium.’ and ’Moon dist.’ are the fractional illumination of 
the Moon, and its angular distance from WASP-6 in degrees, at the midpoint of the transit. 


Date 

Start time 
(UT) 

End time 
(UT) 

Aobs 

Exposure 
time (s) 

Filter 

Airmass 

Moon 

ilium. 

Moon 

dist. 

Aperture 
sizes (px) 

Scatter 

(mmag) 

2009/06/26 



91 

120 

Bessell R 

1.321.05 

0.271 

160.5 

65,90, 110 

1.215 

2009/08/02 

04:18 


175 

90-120 

Bessell R 

1.28 -t 1.44 

0.934 

59.6 

27, 40, 70 

0.939 

2009/08/29 

02:32 


129 

120 

Bessell R 

1.28 -t 1.20 

0.750 

63.8 

28, 40, 60 

0.598 

2010/07/31 

03:51 


193 

80 

Bessell R 

1.45 -t 1.34 

0.686 

42.4 

25, 35, 55 

0.591 


Table 3. Derived photometric parameters from each light curve, plus the interval within which the best fit was searched for using GEMC. 


Parameter 

Search interval 

2009/06/26 

2009/08/02 

2009/08/29 

2010/07/31 

Radius ratio 

0.05 to 0.30 

0.1443 ±0.0055 

0.1444 ± 0.0043 

0.1474 ±0.0017 

0.1454 ±0.0021 

Sum of fractional radii 

0.10 to 0.50 

0.1102 ±0.0060 

0.1109 ± 0.0048 

0.1115 ±0.0025 

0.1114 ±0.0023 

Linear LD coefficient 

0.0 to 1.0 

0.366 ±0.119 

0.397 ±0.116 

0.368 ± 0.077 

0.402 ± 0.067 

Quadratic LD coefficient 

0.0 to 1.0 

0.245 ±0.191 

0.325 ± 0.222 

0.186 ±0.123 

0.192 ±0.134 

Orbital Inclination (degrees) 

70.0 to 90.0 

88.47 ± 0.99 

88.55 ±0.85 

88.33 ± 0.48 

88.36 ± 0.53 

Transit epoch (BJD/TDB) 

±0.5 in phase 

2455009.83622 ± 0.00021 

2455046.80720 ± 0.00015 

2455073.69529 ± 0.00013 

2455409.79541 ± 0.00010 

Longitude of spot (degrees) 

-90 to -1-90 


-26.15 ± 1.52 

21.30 ±0.99 


Co-latitude of Spot (degrees) 

0.0 to 90.0 


78.76 ± 1.58 

72.77 ±1.12 


Spot angular radius (degrees) 

0.0 to 30.0 


12.25 ± 1.40 

12.17 ±0.81 


Spot contrast 

0.0 to 1.0 


0.649 ±0.187 

0.798 ± 0.082 



Table 4. Combined system and spot parameters for WASP-6. The system 
parameters are the weighted means from all four data sets. The spot angular 
size and contrast are the weighted means from the two transits containing a 
starspot anomaly. 


Parameter 

Symbol 

Value 

Radius ratio 

r-p/r* 

0.1463 ±0.0012 

Sum of fractional radii 

Vs + rp 

0.1113 ± 0.0015 

Linear LD coefficient 

ui 

0.386 ± 0.043 

Quadratic LD coefficient 

U2 

0.214 ± 0.077 

Orbital Inclination (degrees) 

i 

88.38 ±0.31 

Spot angular radius (degrees) 

"Tspot 

12.19 ±0.70 

Spot contrast 

pspot 

0.774 ± 0.075 

Stellar rotation period (d) 

Frot 

23.80 ±0.15 

Projected spin orbit alignment (degrees) 

A 

7.2 ±3.7 


epoch and the bracketed quantities represent the uncertainty in the 
final two digits of the preceding number. Fig. 4 and Table 5 show 
the residuals of these times against the ephemeris. The results show 
no evidence for transit timing variations. 

Initially we used the quoted mid-transit time from Gillon et al. 
(2009), but found that this value disagreed with the other 10 mid¬ 
transit times at the 2.2a level. This may be because the value found 
by Gillon et al. (2009) was derived by simultaneously fitting the 
original WASP data plus two incomplete transits from RISE and 
a single complete transit from the FTS. We therefore used the 
same approach as Nikolov et al. (2014) and fitted (using PRISM) 
the archival FTS light curve to determine the mid-transit time. The 
value found using just the FTS data is in better agreement (0.6a) 
with the other ten mid-transit times. Therefore it was decided to use 
the mid-transit time from the FTS light curve in our analysis, not 
just due to the better agreement but also due to the fact that it comes 
directly from a light curve covering a full transit. 


Table 5. Times of minimum light of WASP-6 and their residuals versus the 
ephemeris derived in this work. 

References: (1) Gillon et al. (2009); (2) This work; (3) Dragomir et al. 
(2011); (4) Jordan et al. (2013); (5) Sada et al. (2012); (6) Nikolov et al. 
(2014) 


Time of minimum 
(BJD/TDB - 2400000) 

Cycle 

no. 

Residual 

(BJD) 

Reference 

54425.02167 ±0.00022 

0.0 

-0.00013 

1 

55009.83622 ±0.00021 

174.0 

0.00006 

2 

55046.80720 ±0.00015 

185.0 

0.00001 

2 

55073.69529 ±0.00013 

193.0 

0.00008 

2 

55409.79541 ±0.00010 

293.0 

-0.00000 

2 

55446.76621 ±0.00058 

304.0 

-0.00023 

3 

55473.65438 ±0.00016 

312.0 

-0.00007 

4 

55846.72540 ±0.00045 

423.0 

-0.00028 

5 

56088.71800 ±0.00013 

495.0 

0.00017 

6 

56095.43973 ±0.00017 

497.0 

-0.00011 

6 

56132.41081 ±0.00010 

508.0 

-0.00005 

6 


4.2 Physical properties of the WASP-6 system 

With the photometric properties of WASP-6 measured the physi¬ 
cal characteristics could be determined. The analysis followed the 
method of Southworth (2009), which uses the parameters measured 
from the light curves and spectra, plus tabulated predictions of the¬ 
oretical models. We adopted the values of i, r.plri, and r* + rp 
from Table4, the orbital velocity amplitude TT* = 74.3^\'^ ms“^ 
and eccentricity e = 0.054 ('(q q ]^5 from Gillon et al. (2009), and 
the stellar effective temperature Tes = 5375 ± 65 K and metal 
abundance = —0.15 ± 0.09 from Doyle et al. (2013). 

An initial value of the velocity amplitude of the planet, Kp, 
was used to calculate the physical properties of the system using 
standard formulae and the physical constants listed by Southworth 
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Figure 3. Transit light curves and the best-fitting models of WASP-6. The 
residuals are displayed at the base of the figure. 

(2011). The mass and of the star were then used to obtain the 
expected Ten and radius, by interpolation within a set of tabulated 
predictions from theoretical stellar models. Kp was iteratively re¬ 
fined until the best agreement was found between the observed and 
expected Ten, and the measured r* and expected —. This was per¬ 
formed for ages ranging from the zero-age to the terminal-age main 
sequence, in steps of 0.01 Gyr. The overall best fit was found, yield¬ 
ing estimates of the system parameters and the evolutionary age of 
the star. 

This procedure was performed separately using five different 
sets of stellar theoretical models (see Southworth 2010), and the 
spread of values for each output parameter was used to assign a 
systematic error. Statistical errors were propagated using a pertur¬ 
bation algorithm (see Southworth 2010). 

The final results of this process are in reasonable agreement 
with themselves and with published results for WASP-6. The fi¬ 
nal physical properties are given in Table 6 and incorporate sepa¬ 
rate statistical and systematic errorbars for those parameters which 


Table 6. Physical properties of the WASP-6 system. Where two enorbars 
are given, the first is the statistical uncertainty and the second is the system¬ 
atic uncertainty. 


Parameter 

Value 

At A (Mq) 

0.836 ± 0.063 ±0.024 

Ra (R©) 

0.864 ± 0.024 ±0.008 

logSA (cgs) 

4.487 ± 0.017 ±0.004 

PA (p©) 

1.296 ± 0.053 

Mb ( Mjup) 

0.485 ± 0.027 ±0.009 

-^b ( R^Jup) 

1.230 ± 0.035 ±0.012 

Sb (ms-2) 

7.96 ± 0.30 

Ph ( PJup) 

0.244 ± 0.014 ±0.002 

Tf (K) 

1184 ± 16 

e 

0.0390 ±0.0014 ± 0.0004 

a(AU) 

0.0414 ±0.0010 ± 0.0004 

Age (Gyr) 

Q n +8.0 +4.0 
-12.7 -9.0 


depend on the theoretical models. The final statistical errorbar for 
each parameter is the largest of the individual ones from the solu¬ 
tions using each of the five different stellar models. The systematic 
errorbar is the largest difference between the mean and the individ¬ 
ual values of the parameter from the five solutions. 


5 STARSPOT ANOMALIES 

Two of the light curves, from 2009/08/02 and 2009/08/29, contain 
apparent starspot anomalies (see Fig. 2). Due to a 27 day gap be¬ 
tween the two light curves it is not possible to conclusively demon¬ 
strate that the anomalies are due to the same spot. But if so, the 
stellar rotation period and sky-projected spin orbit alignment can 
be calculated and compared to the values found by Gillon et al. 
(2009), Doyle et al. (2013) and Nikolov et al. (2014). This will al¬ 
low an indirect check on whether the two spot anomalies are due to 
the same starspot. 

Firstly we consider whether the spot could last for a 27 day 
period. On the Sun a spot’s lifetime T is proportional to its size Aq 
following the Gnevyshev-Waldmeier (G-W) Relation (Gnevyshev 
1938; Waldmeier 1955). 

Ao = WT (2) 

where Aq is measured in MSFl (micro-Solar hemispheres) and 
T is in days. Petrovay & van Driel-Gesztelyi (1997) state that 
W = 10.89 ± 0.18 MSH day"^ Kenwood et al. (2010) showed 
that large sunspots also followed the Gnevyshev-Waldmeier rela¬ 
tionship. If the same relationship is applied to starspots then a min¬ 
imum lifetime of 30 days requires a minimum size of 327 MSH, 
or an angular radius of just greater than 1°. Bradshaw & Hartigan 
(2014) argues that the standard solar G-W relation overestimates 
the lifetime of a starspot. Bradshaw & Hartigan (2014) uses tur¬ 
bulent magnetic diffusivity at supergranule size scales to calculate 
the magnetic diffusivity which in turn allows W in the G-W rela¬ 
tion to be re-calculated. Depending on the turbulent scale length 
being used and to have a minimum lifetime of 30 days requires a 
angular radius of 3° to 9° (see fig. 1 Bradshaw & Hartigan 2014). 
The sizes of the starspot anomalies in the WASP-6 light curves are 
greater than 10°, so we conclude that a single spot can last suffi¬ 
ciently long to cause both anomalies, irrespective of the turbulent 
scale length used. 
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Figure 4. Residuals of the available times of mid-transit versus the orbital ephemeris found for WASP-6. The four timings from this work are the cluster of 
three points between the cycle numbers 170-200 and the point close to cycle 290. 


5.1 Starspot anomalies results 

The results from modelling the two spot anomalies suggest that 
they are due to the same spot rotating around the surface of the 
star, as the spot sizes and contrasts are in good agreement and the 
lifetime of a spot this size is much greater than the time interval 
between the two spotted transits. Fig. 5 is a representation of the 
stellar disc, the spot and the transit chord for the two nights of ob¬ 
servations. 

By assuming that the two spot anomalies are indeed caused by 
the same spot, it is straightforward to calculate the sky-projected 
spin orbit alignment of the system. We find a value of A = 7.2° ± 
3.7° from the measured positions of the starspot during the two 
transits. 

It is also possible to calculate the rotational period of the star, 
using the spot positions and an estimate of the number of stellar ro¬ 
tations which occurred between the two transits (see Tregloan-Reed 
et al. 2013; Mancini et al. 2014). Due to the 27 day gap between 
the light curves the star could have rotated N full rotations plus 
47.5° ± 2.5°. If W = 0 then this would imply that WASP -6 has a 
rotation period of approximately 200 days, which is extremely long 
for a main sequence G-star. If A = 1 then the spot has travelled 
407.5° ± 2.5° between the transits, giving a rotational period of 
Prot = 23.80 ± 0.15 d at a co-latitude of 75.8°. This is in excel¬ 
lent agreement with the measurement of Prot = 23.6 ± 0.5 d from 
Nikolov et al. (2014). Combining this with the stellar radius (see 
Table 6 ), the latitudinal rotational velocity of the star was calcu¬ 
lated to be ^( 75 ,go) = 1.78±0.20 km s“^. This is also in agreement 
with V sin I from both Gillon et al. (2009) and Doyle et al. (2013). 
If A = 2 then the spot has travelled 767.5° ± 2.5°, giving a rota¬ 
tional period of Prot = 12.63 ± 0.15 d at a co-latitude of 75.8° (or 
f( 75 . 8 °) = 3.36 ± 0.20kms“^). This agrees with the vsinl from 
Gillon et al. (2009) and Doyle et al. (2013), but not with the Prot 
from Nikolov et al. (2014). The agreement with Gillon et al. (2009) 
and Doyle et al. (2013) is due to the fact that any value of v that is 
found to be greater than v sin / can be considered to agree based 
on the nature of sin I. We conclude that the A = 1 case is much 
more likely than the two alternatives discussed above. 


5.2 Degeneracy of the stellar rotation period 

Whilst there is no clear photometric signal in the SuperWASP 
light curve of WASP-6, Nikolov et al. (2014) were able to mea¬ 
sure a rotation period of Prot = 23.6 ± 0.5 d from photometry 
of higher precision, however, none of the STIS observations de¬ 
tected a starspot anomaly indicating that starspots on WASP-6 A 
are either rare or of low contrast. This is also supported by the 



Figure 5. Representation of the stellar disc, starspot, transit chord and equa¬ 
tor for the two datasets of WASP-6 containing spot anomalies. The axis of 
stellar rotation lies in the plane of the page and in the case of A = 0° points 
upwards. 


upper limit of the photometric variability of about 1% (Nikolov 
et al. 2014). There are also two measurements of usinJ from 
Gillon et al. (2009) (v sin I = 1.4 ± 1.0 km s“^) and Doyle et al. 
(2013) (usinJ = 2.4 ± 0.5kms“^). Both vsinl measurements 
agree with the v found when combining Prot and P* at a co¬ 
latitude of 75.8° to give either ^( 75 ,go) = 1.78 ± 0.20kms“^ 
or W( 75 . 8 °) = 3.36 ± 0.20kms“^. The problem that arises from 
checking measurements of v against vsinl is that due to the 
sin I projection factor any value for v that is found to be greater 
than V sin I can be considered to agree. A second unknown is the 
amount of differential rotation that is experienced by WASP -6 A. In 
the absence of any differential rotation the single full rotation value 
of Prot = 23.80 ± 0.15 d would lead to an equatorial rotational ve¬ 
locity of w = 1.84 ± 0.20 km s“^. This result agrees again with the 
V sin I value from both Gillon et al. (2009) and Doyle et al. (2013). 
Our results from PRISM do show though that the two starspot posi¬ 
tions are only approximately 10° from the stellar equator. As such 
the effect from differential rotation would be small, so any large 
divergence of v from v sin I would imply that I 90°. 

WASP -6 A has T^s = 5375 ± 65 K (Doyle et al. 2013) so 
is a cool star (Teg < 6250 K). The trend seen between host star 
TeflS and projected orbital obliquity (see Fig. 6 ) suggests that the 
orbital rotation axis of WASP -6 b should be aligned with the stellar 
rotation axis of WASP -6 A. For this to be true then I would have 
to be « 90°, and thus sin / « 1. If this is the case then the value 
'^( 75 . 8 °) = 3.36 ± 0.20kms”*^ no longer agrees with the vsinl 
from either Gillon et al. (2009) or Doyle et al. (2013). This supports 
the supposition that the rotation period of WASP -6 A is Prot = 
23.80 ± 0.15 d. Brown (2014) calculated the stellar rotation period 
of WASP -6 A to be Prot = 27. d from Gaussian distribution 
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sampling of vsinZ, i and Ra, which is further evidence for the 
conclusion that the rotation period of WASP-6 A is P^ot = 23.80 ± 
0.15 d. 


6 DISCUSSION AND CONCLUSIONS 

We have determined the physical properties of the WASP-6 plane¬ 
tary system (Table 6) based on four new high-precision transit light 
curves, finding values which are consistent with and more precise 
than those in the literature. We find the mass and radius of the host 
star to be 0.836 ± 0.063 M© and 0.864 ± 0.024 R©, respectively. 
For the planet we find a mass of 0.485 ± 0.027 Mjup, a radius of 
1.230 ±0.035 Rjup and a density of 0.244±0.014 pjup. These re¬ 
sults also serve as a secondary check for the accuracy of the PRISM 
and GEMC codes. By studying the individual results for each of the 
four transits (see Table 3) it can be seen that the system parame¬ 
ters from each light curve agree within their la uncertainties. This 
shows that PRISM can retrieve reliable photometric properties from 
transit light curves containing starspot anomalies. 

The four transits of WASP-6 were modelled using PRISM and 
GEMC. Two of the transits contained a starspot anomaly but are 
separated by 27 days. Whilst it is not possible to prove that the two 
spot anomalies are caused by the same starspot, the available evi¬ 
dence strongly favours this scenario. The results from PRISM show 
that the angular size and contrast of the starspot in both light curves 
agree to within 0.05cr and 0.73 (t, respectively. As with WASP-19 
(see Tregloan-Reed et al. 2013), only part of the starspot(s) is on 
the transit chord (Fig. 5). Because the light curve only holds infor¬ 
mation on what is happening inside the transit chord then a likely 
scenario is that the planet is passing over a band of smaller starspots 
which form an active region on WASP-6. In this active region, there 
could be a number of starspots each with sizes much less than 1° 
and therefore lifetimes shorter than 30 days (see Section 5). Future 
observations may allow changes to be seen in the overall contrast 
from the starspot region. In either case as a whole the region would 
remain a similar size and shape over a 27 day period. 

In the case of a single large starspot, Pspot (Tabled) and R* 
(Table 6) can be combined to find the starspot radius. We find 
Rspot = 127902 ± 11102 km, which equates to approximately 
4.5% of the visible stellar surface. This value is similar to starspots 
found on other G-type stars (Strassmeier 2009). 

If the two starspot anomalies are assumed to be generated by 
the planet crossing the same starspot then it is possible to calculate 
the latitudinal rotation period of WASP-6. It was found that either 
Prot = 23.80 ± 0.15 d or Prot = 12.63 ± 0.15 d at a co-latitude of 
75.8°. These calculations assumed that WASP-6 had made either 
one or two full rotations prior to the difference seen in the light 
curves. 

Even without knowing the number of full rotations that 
WASP-6 completed between the two spotted light curves, if the 
starspot anomalies are due to the same spot then the sky-projected 
spin orbit alignment A of the system can be measured. We find 
A = 7.2° ± 3.7°. This result agrees with, and is more precise than, 
the previous measurement of A using the RM effect (A = 11° 
Gillon et al. 2009). A gives the lower boundary of the true spin-orbit 
angle, tp. As stated by Fabrycky & Winn (2009), finding a small 
value for A can be interpreted in different ways. Either ip lies close 
to A and the system is aligned, or pj lies far from A and the system 
is not aligned. As discussed in Section 5.1 because the spot is close 
to the stellar equator then it could be assumed that the change in v 
at the equator due to differential rotation would be small. Coupled 



Figure 6. |A| against Teff for 83 transiting planets from TEPCat including 
WASP-19 and WASP-6. The green and red datapoints are WASP-6 (left) 
and WASP-19 (right). The green datapoints represent values from the lit¬ 
erature (WASP-6: Gillon et al. 2009; WASP-19: Hellier et al. 2011) and 
the red datapoints represent the values found from this work and Tregloan- 
Reed et al. (2013). The trend in the data suggests that cool host stars harbour 
aligned systems. 


with the uncertainties measured in v sin I from both Gillon et al. 
(2009) and Doyle et al. (2013) it is plausible that sin I « 1 and 
therefore t/i ~ 7° if R-ot = 23.80 ± 0.15 d. As a consequence we 
have two different scenarios: an aligned system with a slowly rotat¬ 
ing star or a misaligned system with a rapidly rotating star. Taking 
into account the Tefr of WASP-6 A and the statistical trend seen in 
misaligned systems it is more probable that the WASP-6 system is 
in fact aligned, suggesting ip ^ 7° and Prot = 23.80 ± 0.15 d. It 
would be desirable to observe consecutive transits of WASP-6 in 
an attempt to definitively identify multiple planetary crossings of a 
single starspot and to precisely determine Prot, A and potentially ip 
of WASP-6. 

If the starspot anomalies are due to the same starspot, A = 
7.2° ± 3.7° and there is no direct evidence for a spin-orbit mis¬ 
alignment in the WASP-6 system. With potentially a low obliquity 
and a cool host star, WASP-6 seems to follow the idea put forward 
by Winn et al. (2010a) that planetary systems with cool stars will 
have a low obliquity. It also lends weight to the idea that WASP-6 b 
formed at a much greater distance from its host star and suffered or¬ 
bital decay through tidal interactions with the protoplanetary disc 
(i.e. either Type I or Type II disc-migration. Ward 1997). 

At present there are 83 transiting planets with published A 
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values^. The A values for WASP-6 (this work) and WASP-19 
(Tregloan-Reed et al. 2013) were updated and a plot of A against 
Teff was created (see Fig. 6). To remove any ambiguity in the plot 
due to negative values of A, we plot its absolute value. It can be 
seen that a large proportion (75 %) of cool stars (Teg < 6250 K) 
are in aligned systems, while the majority (56 %) of hot host stars 
have misaligned systems. This trend supports Winn et al. (2010a) in 
that cool stars with hot Jupiters will have low obliquities. This trend 
can also be explained by the time required for the system to align. 
Hot stars will have thinner convective zones and will therefore take 
longer to align the photosphere with the planetary orbit. Because 
of this, by examining A of hot stars a greater proportion will have 
misaligned systems compared to cool stars where the alignment 
process is much shorter and so will have a higher proportion of 
aligned systems. Cool stars also live longer so the ones that are ob¬ 
served are on average older. They have therefore had more time for 
tidal effects to work (Triaud 2011). 

By determining A and tp of the planetary system it is possible 
to begin to understand the primary process in the dynamical evolu¬ 
tion of the system. The RM effect can be used to ascertain a value 
for A. One limitation of this method though is from an excess RV 
jitter (stellar activity e.g. starspots). Therefore, the use of the RM 
effect either requires magnetically quiet stars or the transit chord 
of the planet to bypass any active latitudes on the stellar disc. The 
opposite is true when using starspot anomalies in light curves to de¬ 
termine A. Due to this the two different methods complement each 
other in probing the dominant process in the dynamical evolution 
of transiting planets. It should be noted that in both the cases of 
WASP-19 and WASP-6 (see Fig. 6) the measured uncertainty in A 
is much smaller than measured using the RM effect. This indicates 
that the starspot method to measure A is superior to the RM effect 
in terms of reduced uncertainty in measuring A. However, as was 
shown in observing WASP-50 (see Tregloan-Reed & Southworth 
2013), the starspot method does not always work in terms of ob¬ 
taining transit light curves affected by a starspot anomaly. The RM 
effect does have a high success rate in measuring a value of A but 
rarely achieves a similar precision. 
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